Frobenius Numbers by Lattice Points

Integer Optimization and a Classic Computational Problem

Daniel Lichtblau, Wolfram Research, Inc.
Stan Wagon, Macalester College, St. Paul, Minnesota
Joint with
     David Einstein (sdcorp) & Adam Strzebonski (WRI)

Abstract

The Frobenius number g(A) of a finite set of positive integers is the least integer that is not representable as a nonnegative combination of entries in A. For example, Chicken McNuggets come is sizes of 6, 9, and 20. One cannot buy 43 McNuggets, but any larger target is representable. Thus g(6, 9, 20) = 43. Finding the Frobenius number, and solving Frobenius instances (finding a solution to X·A=T in nonnegative integers X) are closely related to more general problems in integer optimization. For example, minimizing X·A-T subject to the constraints that this quantity is nonnegative and each entry of X is nonnegative will solve the instance problem, since the target T is representable iff 0 is the minimum.

In this talk we will show how a study of a certain lattice of integer vectors leads to a geometric object — the fundamental domain — that encodes enough information to find the Frobenius number. Getting this information about the domain requires developing several algorithms, many of which have had an impact on some of Mathematica's kernel functions, such as Reduce and Minimize. Letting n be the number of integers in A, we can find g(A) (with rare exceptions) when n ≤ 11, but with no restriction on the size of the entries in A. The algorithm seems to have complexity that, when n is fixed,  is close to quadratic on average. The previous best algorithm for this problem was polynomial-time in the input size, but with polynomial degree (n^n). Here is an example with n = 5 (in version 5.2).

Needs["NumberTheory`"]

A = {10000000001, 20100000000, 30000000017, 40000000039, 
          50999000101};
g = FrobeniusF[A]

FrobeniusInstance[A,  g + 100]

The main points are best illustrated when A = (a, b, c, d). Letting (Subscript[ℤ, a]) be the integers mod a, there is a natural map \(TraditionalForm\`h : ℤ\^3 -> ℤ\_a\) given by h(X)= X.(b, c, d) (mod a). The kernel of h is the lattice L consisting of triples X such that X.(b, c, d) is divisible by a. The classic relationship (ℤ^3)/
          L ≈ (Subscript[ℤ, a]) raises the question of finding a geometric form for (ℤ^3)/L; by this we mean choosing a representative vector for each equivalence class. We do this in a natural way, and call the resulting set of a vectors 𝒟, for fundamental domain; it sets in the first octant and includes the origin. There is then a tiling of (ℤ^3) induced by translating copies of 𝒟 by vectors in L. Here is picture of the domain — the gray cubes — when A = (100, 239, 543, 609).

Code for Art

<< FrobeniusNumber`;
<< Geometry`
tet = Polygon[Vertices[Tetrahedron][[#]]] &  /@  
          Faces[Tetrahedron];

    
A = {100, 239, 543, 609};
InteriorElbow[A];
cor = Corners[A]; elb = Elbows[A];
fd = 
      Block[{$DisplayFunction  =  Identity}, 
        FundamentalDomain[A, ShowCorners -> False, 
          ShowElbows -> False,  
          Colors -> (Red &  /@  Range[4])]]; 
    im1 = Show[fd, 
        Graphics3D[{{RGBColor[0.3, 0.8, 
                0.3], Cuboid [# +  .15 -  .001,  # +  .85 +  .001] & /@  Delete[cor, -2]}; 
{RGBColor[1, 1, 
                0],  Cuboid [# +  .15 -  .001,  # +  .85 +  .001] & /@  elb}; 
{GrayLevel[ .95], RGBColor[1, 1, 0],  
              EdgeForm[
                GrayLevel[
                  0]], (tet /.  {x_?NumericQ,  y_, 
                          z_}  :>  {x, y, 
                              z}  .3  +  #  +   .4)  & /@  
                elb}, 
{RGBColor[ .7, 0.2, 0.2];  
              RGBColor[0.2, 0.4, 
                1], Cuboid [# +  .15 -  .001,  # +  .85 +  .001] &@ 
                cor[[-2]]}; {}}],  
        DisplayFunction  ->  $DisplayFunction];

The trick then is to capture concisely the geometric structure of the domain, so that it can be easily described even when it has (10^100) entries. We do this by getting the elbows of 𝒟, by which we mean vectors E that are not in 𝒟, but such that subtracting 1 in any nonzero coordinate yields a point that is in 𝒟; these are the nine tetrahedra in the diagram. It follows that 𝒟 is obtained by deleting the cones defined by all the elbows. Having the elbows in hand allows us to get the corners of 𝒟, and one of the corners yields the Frobenius number of A.

Our algorithm for doing this is complicated, with many steps. But it appears that, when n is fixed, the average complexity is close to quadratic. Other consequences of our work are a very fast algorithm when n = 3 (it can handle million-digit inputs in under 80 seconds) and some new theoretical results about the Frobenius number of special sequences.

Solving Frobenius Instances

An important part of our algorithm involves solving linear diophantine equations to give important information about the domain 𝒟. Such problems are also of importance in their own right, with applications in computational number theory, operations research, and elsewhere. A particular application, closely related to the problem at hand, is solving a specific Frobenius instance: Given our set A and a target value T, we find a nonnegative X that solves X·A = T. We will discuss how this is done efficiently in Mathematica, and mention how similar diophantine equations relate to finding Frobenius numbers.

The Frobenius Problem

Postage-Stamp Problem

Given positive integers (Subscript[a, 1]), (Subscript[a, 2]), …, (Subscript[a, n]) and a target T find nonnegative integers (Subscript[x, i]) so that \(TraditionalForm\`\(x\_1\) 
            a\_1 + \(\(·\)\(·\)\(·\)\(+x\_n\)\) a\_n = T\) (or prove impossible).  (Same as X·A=T)

Frobenius Number

The number of NONrepresentables is always finite, and so:

The Frobenius number g((Subscript[a, 1]), (Subscript[a, 2]), …, (Subscript[a, n])) is the largest target T not representable as above. I.e., the largest number not in the semigroup over ℕ generated by A.

Example by Hand

g(2, 5) = 3, since one cannot get 3, but it is easy to get 2 + 2 = 4,2 + 3 = 5, 6, 7,…

Chicken McNuggets

A = {6, 9, 20};
X = FrobeniusInstance[A,  44]
X . A

FrobeniusInstance[A,  43]
FrobeniusNumber[A]

Examples

Big Examples: Big n

A = TestData[(10^5), 20]
g = FrobeniusNumber[A]

X = FrobeniusInstance[A, g + 1]
X . A
FrobeniusInstance[A, g]

Big Examples: Big (Subscript[a, i])

A = TestData[(10^20), 6]
g = FrobeniusNumber[A]

X = FrobeniusInstance[A, g + 1]
X . A

FrobeniusInstance[A, g]

Real World Postage Problems

If you have five 39¢ stamps and some 10¢ and 24¢ stamps and need to make $2.34.

{x, y, z}/.FindInstance[
     {x, y, z} . {10, 24, 39} == 234&&
    0 ≤ x&&0 ≤ y&&0 ≤ z ≤ 5,
{x, y, z},Integers]

I wrote USPS that they should add this to usps.com; would be quite useful.

Danny will talk more about finding solutions to such equations.

Our Results

History

g(a, b) = a b  -  a  -  b  (Sylvester, 1880)

g(a, b, c) by a quadratic-time algorithm (Greenberg, Davison, 1990)

g(A) poly-time algorithm known when n is fixed (Kannan), but very complicated — O((((log a))^((n^n)))). No practical value.

From now on we write A = {a, (Subscript[b, 1]), (Subscript[b, 2]), …, (Subscript[b, m])}, so m = n - 1.

Graph methods: Big n, small a (Beihoffer, Nijenhuis, Hendry, Wagon, Elec J Comb, 2004)

Make a graph with a vertices and study shortest paths. There are fast algorithms for this, but it is limited by the size of a since O(a) is needed for time and for space. 

Works for a ≤ (10^7), with essentially no bound on n. Included in Mathematica a couple of years ago and is very fast in its domain.

Lattice methods: Big a, small n (Einstein, Lichtblau, Stzrebonski, Wagon)

Get g(a, b, c) in time that is O((((log a))^(1 + ϵ))) (called softly linear); we can handle million-digit inputs!

Get g(A) in average-time that is O((((log a))^(2 + ϵ)))  (soft quadratic), but performance only conjectured, and constant rises with n, so that n ≤ 11. But a~(10^100) is fine.

Improved Frobenius instance solving; Danny will address.

A = TestData[(10^1000), 3];
Timing[FrobeniusNumber[A]] // Short

A = TestData[(10^100), 4];
Timing[FrobeniusNumber[A]]

Our Lattice Methods

Given A  =  (a, B) where B is an m-vector, there is a natural homomorphism to integers mod-a.

	\(TraditionalForm\`h : ℤ\^m -> ℤ\_a\)    by     h(X)=X·B (mod a)

Kernel of h: vectors X with X·B is divisible by a. This is a lattice L.

Points below form the lattice for A = (33, 89, 147): points (x, y)·(89, 147)≡0 (mod 33)

h is surjective, so (ℤ^m)/L is isomorphic to (Subscript[ℤ, a]). So we can choose representatives in (ℤ^m) for each {0, 1, 2, …, a - 1} — i.e., for each residue r < a find X so that X·B≡r (mod a). Then 𝒟, the set of representatives, has the property that 𝒟 tiles (ℤ^m) when translated by vectors in the lattice L.

This is standard modern algebra. What we have done is made it concrete.

1. For each residue class, we choose the point in 𝒟 so that it is in the first quadrant and has minimal weight (B·X)

2. We relate the Frobenius number to the geometry of 𝒟 (above were known)

3. We look at pictures of 𝒟 and figure out how to characterize it, even if it has (10^100) elements, with just a few vectors.

4. We figure out how to compute those vectors.

The Lattice

Here is a simple case A = (33, 89, 147). The black dots are the lattice points. The domain 𝒟 is the set of 33 squares. The three yellow disks are the "elbows". The amount of blackness at point X is X·(89, 147), the weight of the point.

Tiling

Definition: Elbow: a vector E not in the domain, but such that moving down 1 in any coordinate that is nonzero is in the domain. An axal elbow is an elbow on one of the m axes. 

Definition: Corner: a vector C in the domain but such that adding 1 in any direction goes outside of the domain.

Theorem: Corner C of largest weight yields the Frobenius number by subtracting a.

Finding the axial elbows is an optimization: 9  b is smallest such that 9  b == α a + x 147 with α ≥ 1 & x ≥ 0.

In general, minimize k so that   \(TraditionalForm\`k\ b\_1 = α\ a + \ \(x\_1\) 
            b\_2 + \(\(·\)\(·\)\(·\)\(+x\_m\)\) b\_m\) and α ≥ 1 (gives weight condition) and (Subscript[x, j]) ≥ 0 and k ≥ 1. This is an INTEGER LINEAR PROGRAMMING problem (ILP).

The domain one dimension up.

Here are the elbows and corners (heaviest corner in blue).

Finding all elbows is the crux of what we accomplished. Their number does not, on average, rise as a rises (though it does increase with n). In other words, for fixed n,the shape of the domain does not change much. Indeed, this shows the power of having high precision: understanding such facts helped us formulate the general algorithm. Here is the domain with (10^100) points!!

Domain with (10^100) points

Code

SeedRandom[4];
A = TestData[10^100, 4];
elb = Elbows[A];
cor = N@Corners[A]/10^32;

    
\(\(Show[
        Graphics3D[{EdgeForm[Thickness[ .01]], 
            Red, \(Cuboid[{0, 0, 0}, \ #] &\)\  /@ \ Sort[cor]}], \ 
        ViewPoint\  -> \ {1.2, 1.2, \ 1.2}, \ 
        AxesLabel -> {x, y, z}, \ Lighting -> False, \ 
        Boxed -> False, \ 
        PlotLabel -> 
          StringForm["a=``, elbows and corners: ``, ``", \*""(10^100)"", Length[elb], \ Length[cor]], \ 
        TextStyle -> \(FontSize -> 16\)];\)\)

SeedRandom[4];
A = TestData[10^100, 4];
elb = Elbows[A];
cor = N@Corners[A]/10^32;
lines[Cuboid[{a_, b_, c_},  {d_,  e_, f_}]] := 
      Line  /@  {
{{a, b, c},  {d, b, c}} ,  {{a, b, 
              c},  {a, e, c}}, {{d, e, c},  {a, e, c}}, {{d, e, c},  {d, 
              b, c}}, 
{{a, b, f},  {d, b, f}} ,  {{a, 
              b, f},  {a, e, f}}, {{d, e, f},  {a, e, f}}, {{d, e, 
              f},  {d, b, f}}, 
{{a, b, c},  {a, b, 
              f}} ,  {{a, e, c},  {a, e, f}}, {{d, e, c},  {d, e, 
              f}}, {{d, b, c},  {d, b, 
              f}}
 }
\(\(Show[
        Graphics3D[{EdgeForm[], \({Hue[Random[]], Red, 
                  cubes = Cuboid[{0, 0, 0}, \ #\  .95]} &\)\  /@ \ 
              cor, 
cubes = \(Cuboid[{0, 0, 
                      0}, \ #\  .95] &\)\  /@ \ 
                cor; 
{Black, Thickness[0.007], \ 
              cubes\  /. \ 
                Cuboid\ [a_, \ b_]\  :> \ 
                  lines[Cuboid[a\ \ , b\ \ ]]}
}], \ 
        ViewPoint\  -> \ {1.2, 1.2, \ 1.2}, \ 
        AxesLabel -> {x, y, z}, \ Lighting -> False, \ 
        Boxed -> False, \ 
        PlotLabel -> 
          StringForm["a=``, elbows and corners: ``, ``", \*""(10^100)"", Length[elb], \ Length[cor]], \ 
        TextStyle -> {FontFamily\  -> \ "Times", 
            FontSize -> 16}];\)\)

Once the elbows are in hand, we can use an elbow-to-corner algorithm (worked out by Einstein) to get the corners, and then the largest-weight corner.

So how to find the elbows? Every elbow comes from a protoelbow. Look at (9, 0) in example (33, 89, 147). It is not in the domain 𝒟. (0, -5) is the point in 𝒟 in the same congruence class. So (9, -5) is the vector that tells the whole story. It is always like this. And (9, -5) is a lattice vector: (9, -5)·(89, 147)=66 ≡ 0. So we need only find all lattice vectors in the bounding box — the green points. Adam Strzebonski figured out how to do this by a simple recursive algorithm (Reduce, which is how I first realized it could be done).

Example

SeedRandom[1]; A = TestData[(10^6),  5]
AxialElbows[A]
ProtoElbows[A] // Short
PreElbows[A] // Short
Elbows[A]
Sort[
      Corners[A],  # . Rest[A]  ≤  #2 . Rest[A] &]

{-67, -36, 12, 23} . Rest[A]

So (0, 0, 12, 23) is a potential elbow.

{86, 6, 23, 19} . Rest[A] - (10^6)

FrobeniusNumber[A]

This all works up to n = 11, but there are 500,000,000 protoelbows, so it takes a while.
There are many other tricks we use to speed things up, but this is the basic idea.

The Bad News

In rare cases (not in random inputs), the number of elbows is so large that our method is doomed. Szekely and Wormald have a formula: (2  (w^2), 2  (w^2) + 1, 2  (w^2) + 2  w, 
        2  (w^2) + 2  w + 1); domain as shown below.

So for a = (10^100) there can be (10^50) elbows. This is why we do not claim a complete solution to the Frobenius problem when n = 4. Yet there is so much structure in these examples that perhaps there is a way.

Efficiencies

There are many tricks used to speed things up in practice. My favorite:

It turns out there is always one (or zero) elbows strictly in the positive orthant. Called minimal elbow. Its weight is the minimal weight. If we consider only the protoelbows having weight less than this minimal weight, that is good enough to get the a set of elbows that might be "rich enough". Any set of elbows can be used to get the corresponding set of corners, and then we can look at the farthest one. It might not be the real Frobenius corner, but its weight (less a) gives an upper bound on g(A).

But whenever we have an upper bound g(A) ≤ U and we suspect g(A) = U, we can certify it as follows:

We know every target past U is representable because g(A) ≤ U.
We can check X·A==U by instance-solving that Danny Lichtblau will talk about. If there is no X, this proves g(A) = U.

SeedRandom[1]; A = TestData[(10^15), 7]

Looking at only some of the elbows yields the upper bound 7560169718000442805. The following proves it is g (A).

FrobeniusInstance[A, 7560169718000442805] // Timing

Bounds

When n ≥ 12, we can get bounds quite quickly. These bounds are much better than what is in the literature.

Here is data for n = 10 and a = (10^10); numbers are base-10-log. In this case we know g(A) = 3.5  …×(10^12). Bounds put it between 3.16 ×(10^12) and 7.76×(10^12).

Typical method: Use the axial elbows only. The corner they define is an upper-bound, and they give us points in the domain, which give lower bounds on the corner.

({{Lower Bound Method, Lower Bound, Upper Bound, Upper Bound Method}, {One random corner, 12.50, 12.89, Initial(7)}, {100 random domain points, 12.49, 13.59, Centerline}, {Axial, 12.43, 13.33, Axial}, {Axial plus minimal, 12.43, 13.30, Axial and minimal}, {Axial,  minimal, and dual axial, 12.43,  ,  }, {Nonrepresentable, 12.30,  ,  }, {Actual value of (Subscript[log, 10])g(A), 12.54,  ,  }})

Application: A New Computer-Derived-&-Proved Theorem

Formulas for the Frobenius number of arithmetic progressions are well known. What about quadratic? I looked at g(a, a + 1, a + 4, a + 9, a + 16) as an example. There were enough patterns that a formula became evident, and then our methods actually led to proof!! Look at one mod-9 congruence class only (a = 9  k); others are similar.

k = (10^6); 
    A = {9  k, 9  k + 1, 9  k + 4, 9  k + 9};
FrobeniusNumber[A]

k = (10^20); 
    A = {9  k, 9  k + 1, 9  k + 4, 9  k + 9};
FrobeniusNumber[A]

k = (10^6); A = {9  k, 9  k + 1, 9  k + 4, 9  k + 9}; 
    ProtoElbows[A]
Elbows[A]
Corners[A]
FrobeniusNumber[A]

k = (10^20); 
    FrobeniusNumber[{9  k, 9  k + 1, 9  k + 4, 9  k + 9}]

Clear pattern: axial elbows are (4,  5, (10^6)); there are three other elbows; corners are simple. 
A proof must start with the protoelbows. Need to prove that (-2, 5,  (10^6) - 2) or in general (-2, 5,  9  k - 2) is a lattice point of positive weight.

Clear[
      k];  {-2, 5, 9 k - 2} .  {9  k + 1,  9  k + 4,  9  k + 9} // 
      Simplify

Aha: Positive and divisible by 9  k. QED. Similar for others. This proves that the elbows might be
 (4, 0, 0)         (0, 
          5, 0)      (0, 0, k)      (1, 2, 
          0)      (3, 0, 1)        (0, 3, 
          k - 1) 

To prove they are exactly these we need to show that the cones determined by these form a region of the right volume (9  k), for then there can be no other bites taken out. 

But the volume is a linear function of k. So ince we know it is correct in two cases, it is proved correct in all cases. It follows that the elbow set is complete and this gets us a proof that the corner set is as we think:

(3, 1, 0)        (0, 2, 
          k - 1)      (2, 1, k - 1)       (0, 
          4,  k - 2)    

with the last one being the heaviest.

Theorem.   g(9  k,  9  k + 1,  9  k + 4,  9  k + 9,  
          9  k + 16) = 9  (k^2) + 18  k - 2

Similar in the other 8 congruence classes. 


Lots of other patterns here; very exciting to discover the formulas, even more exciting to get a complete proof. 

Works up to (a, a + (1^2), a + (2^2), …, 
        a + (7^2)) but then memory and time issues hurt the computational approach, though the formulas no doubt exist and can be discovered; getting proofs is where the memory is required. There were some surprises.

Conclusion

A multi-pronged approach using:

graph theory, lattice theory and algebra, 3-D graphics, high-precision arithmetic, linear programming, integer linear programming, and several co-workers with different areas of expertise, has led to 

• algorithms for the Frobenius problem that go far beyond what was known;

• new code for the kernel that allows Mathematica to make these algorithms available;

• new theorems, with proofs by computer;

• new ideas and algorithms in integer linear programming, which Danny Lichtblau will now speak about.

Init Code

<< FrobeniusNumber`
<< NumberTheory`
FrobeniusNumber[A_]  := 
        If[Length[A] ≥ 4, FrobeniusF[A],  
          FrobeniusFLLLMult[A]];